Wavefunction preparation and resampling using a quantum computer 



Alexei Kitaev and William A. Webb 
California Institute of Technology, Pasadena, CA, 91125-3600, USA 



We present an algorithm that prepares multidimensional Gaussian wavefunctions on qubit arrays 
and an application of such wavefunctions to multidimensional resampling, a technique useful in 
quantum digital simulation. 



INTRODUCTION 



Shor's algorithms for factoring and discrete logarithms [T], as well as Grover's search algorithm |2] promise substan- 
tial quantum speed-up for certain important problems. Another exciting application area for quantum computers is 
simulation of physical systems, including molecules, strongly correlated spin systems, nuclei, and QCD. The majority 
of these problems admit a natural description using continuous variables, such as particle coordinates or field modes. 
("Mode" may refer to the field strength at a grid point or to a Fourier component, depending on which representation 
is chosen.) However, we want to perform the simulation using a digital quantum computer, i.e., one composed of 
qubits. In this paper, we address some general issues pertaining to the discretization of continuous variables and 
round-off errors in the quantum case. It is worth noting that continuous variables play an important role in algebraic 
number theory (for example, consider the use of logarithmic coordinates in the proof of Dirichlet's unit theorem [3' .) 
Thus, we anticipate that discretization techniques might also be helpful in constructing new quantum algorithms for 
number-theoretic problems. 

To see why the discretization is not entirely trivial, consider a quantum particle characterized by a single coordinate 
X. To represent it in a quantum computer, we approximate x by a; = j5, where j is a (binary-encoded) integer and 
5 is some small constant (the discretization grid spacing). To simplify the notation, we may rescale x so that 5 — \; 
thus, we approximate real numbers by integers. Now, suppose we want to change this standard grid to a different 
one, which is equivalent to stretching or squeezing the wavefunction and putting it back on the same grid. We call 
this transformation resampling. If x is mapped to /(x), a continuous wavefunction would transform as follows: 



In the simplest case, f{x) = ^, so that x ^ ^ and '4'{x) i— > y/aip{ay). Naively, this appears to be implementable as a 
permutation of basis vectors. However, the discrete version of /, 



is not a permutation! For a < 1 the transformation is a stretching; it maps integers to a subset of integers, leaving 
some gaps. For a > 1 it is a squeezing, therefore two distinct integers may be mapped to the same integer. 

Fortunately, transformation ^ is to be applied to a special class of wavefunction, namely, ones that vary at distances 
much larger than the grid spacing. In this case, a few lowest bits of x are almost in the uniform superposition and 
thus can be easily removed or added again while maintaining quantum coherence. But adding or removing k bits 
corresponds to stretching (resp. squeezing ) by factor 2*^. It is easy to generalize this trick to any scaling factor a, or 
to an arbitrary function /. 

Multidimensional resampling is somewhat more challenging though. Our solution for this case uses multidimensional 
Gaussian states. Of course, the preparation of such states begins with the one-dimensional case. The transition to 
higher dimensions requires some care. One may be tempted to say that a general Gaussian is just a product of 
several one-dimensional Gaussians. However, the implementation of this approach involves the diagonalization of the 
covariance matrix and a linear change of variables, which in turn requires resampling. To avoid this problem, we use 
a special linear transformation, namely shearing, and construct its discrete version explicitly. In the rest of the paper, 
we simply expand this outline. 
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i:{x) i-> ^(y) = (fix)) ^/j{y), where x = f^^{y). 



(1) 



m - L/o-)j = Lf J 



(2) 



PREPARING A GAUSSIAN IN ONE DIMENSION 



The states we wish to prepare are discrete approximations to continuous Gaussian functions. In this section we 
examine these functions and the approximate forms of interest. In one dimension, the Gaussian is just the normal 



2 



distribution familiar from statistics: G{x,<7,ii) = e "^"'^ (up to a normalization factor). Our algorithm is similar 
to the uniform superposition preparation technique and its generalization to efficiently integrable distributions [S]. 
However, the latter procedures begin with dissecting the distribution at xq — 2^~^ and setting the highest bit of x 
according to |V'(a;)P dx. In comparison, we first set the lowest bit and proceed with generating a discrete Gaussian 
on either even integers or odd integers. This method is slightly more efficient since it does not involve the Gauss error 
function, but rather, the Jacobi theta function, which is easier to calculate. 
If we treat the Gaussian wave function as continuous, it has form 

V'a,p(a;) = — e 2.2 (3) 

so that 

= ^e'^ (4) 

and 

> 

If we consider it to be a function of an integer value, we have 

= e~'2^' (6) 

where the function / is related to the third Jacobi theta function {63) and is defined as 

/(a,M)= £ (7) 

n— — 00 

SO that 

00 

2 — — 00 

If a ^ 1, then /(cr, fi) is approximately equal to e ^ dx = (jy^. We now define the function 

00 

E ^'i^+j-^') (9) 

Note that this function is periodic with period 2''. This allows us to define the quantum state 

2"-l 

\Uf^,N) = E N) (10) 

1=0 

For ^,2''^ — fj, ^ a ^ 1, this state is very close to a Gaussian with mean /i and standard deviation a. A recursive 
description of the state, 

l^<T,M,w) = Iff ,f ^Ar_i) (8)cosa|0) + ,^^^_^) (K)sina |1) , (11) 

where 



a = cos 



- (5. f )//('-"). (12) 



leads almost immediately to an algorithm for preparing it. Note that the sum represented by /(cr, /x) converges, so 
we can approximate it with a finite number of terms and keep some number of digits of the answer. However, the 
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convergence is poor when a is large. When a is small (which is the case at a late stage of our recursive algorithm), 
a few terms of the sum will dominate, and only these need be calculated and stored. For large a we use f{(j,n) = 

ay^exp (^—^^ f —i-;^) to avoid having to sum an exponentially large number of terms. 

The algorithm recursively calls the following subroutine, using the input data (ctojMo) = ^^id some number 

of supplementary qubits. Note that the N qubits are indexed to iV — 1, where the 0-th qubit corresponds to the 
rightmost qubit in a binary representation. 



(a) Calculate and store a (see eqn. (12l ) and a/27r in a k-bit binary representation on supplementary qubits. 

e 0-th qubit 

cos a — sin a 

R{a) = (13) 



(b) Apply the rotation operator R(a) to the 0-th qubit 

cos a — sin a 



sm a cos a 

using the bits of a/27r as described below, 
(c) Uncompute a and a/27r. 

2. Compute (f7i,/ii), where iTi=cro/2 and /ii — ij.o/2 if the previously rotated qubit is |0) and fii — (fiQ — l)/2 if 
it is Note that all pairs (a^, /ii) are stored in supplementary memory until a cleanup phase, see below. 

3. On the remaining — 1 qubits prepare state l^o-i./^i.Af-i)- 

Note that for the last ((N-l)th) qubit, we proceed only through step 1 (c), as after this qubit is rotated we do not 
need another pair of parameters. After this last qubit is rotated, right before termination, a last series of steps is 
performed in which all stored pairs ((7^, /i^) are cleared. To do this, the machine erases (cttv-i, pln-i) by recalculating it 
using {(Jn-2, aud the last output qubit, performing a second controlled XOR to the region where (ctat^i, fiN~i) 

is stored, thus clearing (crjv-i, M-ZV-i) from memory. The process is repeated on (crjv-2, Ai7V-2), (c7V-3, Mw-3) and so 
on until the supplementary memory is cleared. 

Since quantum recursion can be counterintuitive and difficult to describe verbally, this algorithm is perhaps easier 
to understand when presented in circuit form. The figure shows the gist of the algorithm as well as some of the 
workings of the supplementary qubits. Note that the double line across the top denotes the classical part of the 
computer, which in this scheme functions only to transfer the specifications of the Gaussian to the quantum part at 
the end of the process to clear memory. The single lines with forward slashes (/) overlying their left ends represent 
multiple qubits in the supplementary memory. Note that each such line represents the same array at different stages 
of the computation. We deviate slightly from conventional notation by introducing the downward arrow conditional, 
which represents a generalized form of controlled operation. 

The downward arrow terminating on the XOR symbol on the supplementary array directly above the ith qubit 
denotes that the previously used pair (ct;-!, fJ-i-i) and the (i— l)th qubit (which was just rotated) are used to control the 
calculation of {ai,fj,i) and a^. Note that classically, ai is a rotation angle equal to a/ /(cri_i/2, /ii_i/2)//(cri_i, /Ui-i) 
if the {i — l)-th qubit is |0) and \//(CTi_i/2, (^i_i — l)/2)//(CTi_i, /Xi_i) if the (z — l)-th qubit is |1). However, since 
the previous qubit is always in a superposition state, ai is always entangled with the previous qubits. 

The downward arrow terminating on the i-th rotation operator box indicates that the rotation angle is specified 
by the controlling qubits, specifically, the qubits in the supplementary containing ai (which is shown directly to the 
right of the downward arrow). This rotation is actually a series of standardized rotations, each controlled by a single 
qubit of the register storing a^, as described below. 

Due to finite computational resources, the rotation operator must be implemented approximately. Note that since 
a is a rotation angle between and 2tt, so < a/27r < 1 and we can approximate in fc-digits of binary notation that 
a/27r « J2i=i ft so that a « J2i=i '"'jSt- Note that the ai are either 1 or 0, as this is a binary representation. The 
rotation operator is realized as a sequence of at most k standard rotations as R{a) « i?(7r/2'^~^)"'= • • • R{tt/2)°-^R{tt)'^^ . 
The precision of this operation on the entire array (assuming the standard rotations are exact) is thus S — 0{N2^^) 
and the size of the circuit required to execute the algorithm is polynomial in n + log(l/(5). 



PREPARING A GAUSSIAN IN MULTIPLE DIMENSIONS 

In S dimensions, the general Gaussian function is G{x) = Ce~^ Ax/2 ^ where C is a constant, A is a real symmetric 
positive definite matrix {A^=A, and Ax > Vx 7^ 0), and x is vector with S components, x ~ x' — fl where x' 
denotes position in S space and /2 is a vector denoting the mean of the distribution. 
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FIG. 1: Schematic circuit for preparing a Gaussian wavefunction. The top part of the figure shows the general scheme of the 
recursion. The lower figure is a more detailed view of what occurs in a representative block of the process like the boxed section 
of the circuit above. The top and bottom multi-qubit wires represent different regions of the semipermanent supplementary 
memory where all pairs {at, in) are stored until the end of computation. The multi-qubit wires feeding into the black box G 
are the "working" supplementary memory that is inverted and cleared in each step. 
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We assume that the register (call it R) contains Sk qubits, where k is some integer. Thus, we can break the array 
into S subarrays (index them i?j for 1 < z < S) each containing k qubits, and write R=RiR2...Rs, with multiplication 
representing array concatenation. The state space corresponding to R is of dimension 2^'^. The wavefunction we wish 
to prepare has form 

\i>) = Cj2 e-""^"/^ |n) (14) 

In this expressions = 7r~'^/^\/det A and n G , B = [—2*^"^, 2*^"^ — 1], where negative integers are represented 
using two's complement. In the quantum computer, |n) is the basis vector formed by concatenating all binary 
representations of the entries of n. For example, if n = 9, in binary representation, n — 1001, so if we assume k—2, 
n = (10,01)"'" = (2, 1)-^. For the basis state |m) corresponds to the integer m along the ith coordinate axis in S 
space. Thus, each subarray generates one coordinate axis of the space in question. Note that this space is not a state 
space in the system, rather, each integer point in the cube described corresponds to one dimension/basis vector of 
the system's state space. The reason for the distinction will be clear momentarily. 

Our method of preparing the general S'-dimensional Gaussian relies on a special case of the problem. Note that if 
the matrix A in ^^/^ = Ce~^ "^^/^ is diagonal (in which case, we call it D), then G{x) factors into 1-dimensional 
Gaussians as G(x) = YiiLi <^-xp(— rfi.7:|/2) where rf,; is the «th diagonal element of D and Xi the i-th entry of x. 

Thus in this special case, the state to prepare becomes (with n=(ni, ns), where the rij are /c-digit binary integers) 



neB^ fieBSi=l 1=1 VrteB / 

This state is just a product of 1-dimensional Gaussians, which we have already shown how to prepare. We can 
prepare this state by preparing N simple Gaussians on k qubit subarrays. The overall normalization constant (here 
denoted by C) is accounted for by the normalization constants of the Gaussians on subarrays. 

Some simple properties of linear algebra and a change of variables allow us to transform this easily prepared state into 
any arbitrary Gaussian. Any real symmetric positive definite matrix A can be decomposed as ^ = {M'^)^^ D{M^^) 
for some D and M, where D is diagonal and M is upper triangular with I's on the diagonal (in other words, a 
multidimensional shearing matrix). Thus, if we apply transformation n — > ^ = Mn ft = M~^x, then 

exp(-n^i?n/2) = exp ^ ^— = exp — (16) 



The matrix M can be decomposed as the product of shearing matrices with exactly one non-zero off-diagonal 

entry each. Note that M acts numerically on the basis vectors in a classically reversible fashion (provided relevant 
parameters are stored). It is not a unitary operator and is not applied directly to the state space of the system. 
Rather, it maps a vector n to a vector x (and the corresponding basis vector |n) to |a;)) as in the following example. 

Let M be a square, 2 b^ 2 upper triangular matrix with I's the diagonal and a real number a in the upper right 
entry and let n = (ni,n2) • Then, we see that 

M=([ !V ("']=Mn=[ „ 11 ' 1 = 1 - " 1-1 ■ ^ 1 (17) 







(18) 

Thus, M can act on basis vectors numerically by matrix multiplication with rounding, which, provided the speci- 
fications for M are stored, is an efficiently reversible classical operation that can thus be carried out on a quantum 
computer. Applying M to a superposition of states results in a changes of variable of the kind described above, which 
turns our set of independent Gaussians into the general Gaussian, 

C J2 e-"^^"/2 \n)-^C J2 e-^^^^"/2 |^ _ (19) 

HeB^ SeB^ 



Thus, we have an algorithm for preparing the general Gaussian function in N dimensions. 
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RESAMPLING 



One application of the Gaussian states described above is multidimensional resampling. This is more clear after a 
description of resampling in one dimension. The goal is to stretch or contract a wavefunction with respect to some 
function /, see eqn. ([T|. (We assume that / is smooth, monotone, and its derivative is bounded from above and from 
below.) Such deformations are easily accomplished, provided that ip is a, slowly varying function (ip^x) and ip{x ± 1) 
are very close for all x in the range of interest). Assume we have two qubit ensembles, A and B, each containing 
N qubits. The basis states of ensemble A correspond to the x axis, those of ensemble B to the y axis. The desired 
original wavefunction -0 is prepared on ensemble A, while ensemble B is prepared in state |0^). The state of the 
combined system is thus 



Iryo)- ^ V(2;)|x,0^) (20) 

a;=0 

If we want to perform a resampling with respect to f{x) = [|J , then we prepare the uniform superposition state 

TL— 1 

7^ E \J) (21) 



on ensemble B using an algorithm described in Problem 9.4 of [3]. We describe how parameter n is selected shortly. 
After this operation, the whole ensemble is in state 



2™-l n-1 

= ^ ^{X) ^ \X,J) (22) 

' Zn — 



a:— j — — n 

A simple reversible computation can efficiently transform this state into 



2"-l n-1 

)-:7^E^(-)E h ^ +^-> (23) 



V2. , 

Recall that we wish to prepare the state ^{y) — y/ail^ {f{ay)) on an A^-qubit ensemble. Note that if we started with 
that wavefunction prepared on ensemble B with ensemble A initialized to |0^), we could carry out a transformation 
analogous to that above to yield the state 

^ 2"-l an-l 

1^2) = ^ E ^("y) E lL«yJ +j^y)- (24) 



2n 

y— U J — —an 

If we think of A and B as generating coordinate axes x and y, then wavefunctions can be visualized as sets of points 
in the xy plane, specifically, as points occupying the vertices of the infinite grid with integer spacing. In this scheme, 
I772) and 17)2) appear as bands of width 2n about the line y — ^- The amplitudes of basis vectors in \ri2) and |7)2) agree 
on that line; the only difference is that the amplitudes in I772) are constant on vertical strips whereas the amplitudes 
in 1 7)2) are constant on horizontal strips. If n is much larger than the grid spacing but small enough relative to the 
grid length at which ip and ^ vary, the amplitudes in intersecting horizontal and vertical strips are close. Thus we 
have the approximate equality |r/2) = \V2)- This means that by choosing N and n appropriately and preparing the 
state I772), it is virtually equivalent to having prepared the state \ fi2)- If we then apply (to the combined ensemble) 
the inverse of the computation we would have used to prepare \ fi2)i it is as if we prepared and then uncomputed |ry2), 
leaving the combined ensemble in (or rather, arbitrarily close to) the state such that the wavefunction on B is ^{y) 
and the wavefunction on A is jO''^). Thus if we can prepare i^ix) satisfying the necessary smoothness conditions, we 
can use resampling with supplementary qubits to transform the state to one arbitrarily close to ^(j/). 

This particular transformation takes advantage of the fact that uniform superpositions on an interval are invariant 
under stretching or contracting transformation. Multidimensional resampling requires a family of states that is 
similarly invariant under a general linear transformation. Multidimensional Gaussian states satisfy this condition and 
thus may be used for multidimensional resampling. 
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A 2^-1 

FIG. 2: Conceptual representation of states {7)2) and \ri2). The grey region (and the dotted region in the inset) represents the 
basis states of the combined ensemble AB that have non-zero amplitudes. In state \ri2) the amplitudes on vertical strips are 
equal, in state \ f]2) amplitudes arc constant on horizontal strips. If n is much larger than the spacing between basis vectors but 
much smaller than the size of the grid, the amplitudes of basis vectors in a vertical strip will be roughly equal to those of the 
basis vectors of any horizontal strip that intersects it. When this is the case, \fi2) = \r]2) and resampling will work. 



We illustrate with a simple example. Note that in the procedure above, if we prepare on ensemble B the state 

\£.<7,ij.,n) with parameters a and n such that 

\U,,n) ^ E e""^ 1^) = l^-./') = E e""^ 1^) (25) 

i=0 i=—oo 

with /i ^ CT ^ 1 and a = 2n instead of the uniform superposition state \Kn), the procedure still works. Intuitively, 
we prepare a series of vertical, parallel Gaussian-weighted strips centered somewhere above the x axis and shift them 
up along the line y = ^ so that each strip is centered at the line. The result is that the entire grid has non-zero 
amplitudes, but all amplitudes except those immediately about y = f are very, very small, so the result looks very 
similar to Figure 2 and is amenable to resampling as described. 
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